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Abstract 

We  describe  calculations  of  a  plane-wave  pulse  train  of  finite  duration  in  a  linear  lossy 
dispersive  dielectric  half-space  by  means  of  Fourier  integral  transforms.  The  polarization 
response  of  Debye  and  Lorentz  media  to  a  finite  number  of  sinewave  pulses  and  the 
response  of  a  Lorentz  medium  to  a  single  Gaussian  modulated  sinewave  pulse  are  given. 
These  solutions  are  formulated  for  the  case  of  normal  incidence  to  the  planar  boundary  of 
a  dielectric  half-space.  Our  results  show  that  the  Fourier  transform  accurately  reproduces 
precursor  phenomena  previously  observed  for  infinitely  periodic  pulse  trains  using  the 
Fourier  series  method  while  extending  the  analysis  to  aperiodic  waveforms.  These  results 
can  be  computed  to  arbitrary  accuracy  because  the  Fourier  integral  representation  of  the 
propagated  pulse  is  exact,  i.e.  no  analytical  or  physical  approximations  are  used  to  arrive 
at  the  time-domain  solution.  We  believe  that  this  numerical  tool  will  be  useful  in 
developing  physical  intuition  about  the  dynamics  of  pulse  propagation  in  bio-tissue  which 
is  now  made  possible  by  the  availability  of  high-speed  computer  hardware. 

Introduction 

The  propagation  of  wideband  electromagnetic  pulses  in  biological  tissues  has  important 
health  and  safety  ramifications.  A  theoretical  understanding  of  the  dynamics  of  transient 
pulse  propagation  in  dispersive  biological  materials  is  essential  to  the  development  of 
exposure  guidelines  for  non-ionizing  electromagnetic  radiation.  Our  goal  is  to  develop  a 
computational  method  for  wideband  pulse  propagation.  The  method  we  use  is  based  on 
the  Fourier  transform  that  allows  us  to  investigate  the  dynamical  evolution  of 
electromagnetic  pulses  in  linear,  homogeneous,  dispersive  media  by  conducting 
numerical  experiments.  Previous  investigators  have  used  a  Fourier  series  calculation  to 
successfully  predict  the  occurrence  of  so-called  Brillouin  precursors  in  a  water  half-space 
for  incident  pulses  having  sufficiently  short  rise  times  [1].  This  report  documents  the 
computation  of  propagated  waveforms  in  linear  dispersive  half-spaces  characterized  by 
Debye  and  Lorentz  permittivity  models.  These  models  were  chosen  to  demonstrate  the 
Fourier  transform  method  of  calculating  the  transient  response  of  bio-tissues  to  both 
microwave  and  optical  pulsed  signals.  The  essential  difference  between  the  Fourier  series 
and  the  Fourier  transform  is  that  the  spectrum  of  the  latter  is  continuous  and  hence  the 
synthesis  of  the  aperiodic  time-domain  signal  from  its  spectrum  is  accomplished  by 
means  of  integration  instead  of  summation.  The  specific  aperiodic  waveforms  addressed 
in  this  report  are  square-wave  modulated  sinusoids  of  finite  duration  and  the  single 
Gaussian  modulated  sinewave  pulse.  Analytical  and  numerical  aspects  of  the 
computations  are  discussed. 

Formulation 

When  the  frequency  response  of  the  medium  is  available  from  knowledge  of  its  transfer 
function,  then  pulse  propagation  is  straightforward  using  the  Fourier  transform  method. 
We  decompose  the  transmitted  pulse  into  harmonic  waves,  because  harmonic  waves  have 
the  unique  property  of  retaining  their  shape  even  while  propagating  through  highly 
dispersive  media.  However,  the  amplitude  of  the  harmonic  wave  may  decrease  due  to 
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absorption  by  the  medium.  It  is  only  because  harmonic  waves  do  not  change  shape  that 
the  wave  speed  is  precisely  defined.  For  this  reason,  our  discussion  of  waveforms 
propagating  in  dispersive  media  will  be  carried  out  by  resolving  the  waveform  into  a 
superposition  of  harmonic  waves. 

In  a  dispersive  biological  medium,  each  sinusoid  of  frequency  a ,  moves  at  its  own  speed 
and  has  its  own  characteristic  absorption  constant.  Therefore,  our  strategy  for  pulse 
propagation  will  be  to  resolve  the  incident  pulse  shape  into  a  superposition  of  harmonic 
waves.  Propagate  each  wave  of  frequency  0)  for  a  fixed  amount  of  time  at  its  own  phase 
velocity  and  then  resynthesize  the  transmitted  pulse  by  adding  up  all  component  harmonic 
waves  by  means  of  the  inverse  Fourier  integral  while  accounting  for  the  attenuation  of 
each  component  sinusoid  over  the  distance  traveled. 

Frequency  Spectrum 

We  begin  with  a  review  of  fundamental  theory.  Let  the  time-history  of  a  pulse  be 
designated  f(t).  Its  Fourier  spectral  density  F{(0)  is  given  by 

F((0)=  \f(t)e-“°‘dt  (1) 

— OO 

The  Dirichlet  conditions1  are  sufficient  to  ensure  the  existence  of  this  transform  for  all 
physically  realizable  waveforms,  f(t).  The  Fourier  spectrum  F(co )  is  a  complex  number 
used  to  describe  the  amplitude  and  phase  of  a  sinusoid  at  angular  frequency  0) .  The 
complex  amplitude  of  the  Fourier  component  at  angular  frequency  a  is  F{co)do ;  the  area 
under  the  Fourier  spectrum  curve  in  the  interval  do .  This  area  has  the  same  units  as  the 
time-history  of  the  pulsed  electric  field,  i.e.  Volts/m.  Therefore,  the  dimensional  units  of 
the  Fourier  transform  are  different  from  those  of  the  signal  waveform  which  is  being 
transformed.  Specifically,  they  are  the  units  of  the  signal  multiplied  by  time,  i.e. 
(Volts/m)/Hz  or  Volts/m-sec.  The  quantity  F(o)elCCtdo  is  the  sinusoidal  variation  at 
angular  frequency  0) .  The  time-history  of  the  pulse  waveform  can  be  recovered  by 
integrating  its  Fourier  spectral  density  function  or,  in  other  words,  evaluating  the  inverse 
Fourier  transform 

f(t)=  —  f  'F(0))eimd(0 .  (2) 

2  K 

Eqs.  1  &  2  define  the  direct  and  inverse  Fourier  transforms  respectively. 


1  1.  The  waveform  ,/ff)  has  a  finite  number  of  finite  discontinuities, 
has  a  finite  number  of  maxima  and  minima,  and  is  absolutely  integrable,  i.e. 
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Transfer  Function 

The  transfer  function  of  a  linear  time-invariant  system  is  defined  as  the  Fourier  transform 
of  the  impulse  response  of  the  linear  system.  The  transfer  function  of  a  linear  dispersive 
medium  is  the  steady-state  response  per  unit  sinusoidal  input  as  a  function  of  angular 
frequency  co ,  and  depth,  z.  Let  £,(z,tw)  represent  a  frequency  domain  x-polarized  plane 
wave  propagating  in  the  z-direction.  It  is  well-known  that  Ex(z,oo)  must  satisfy  the 
harmonic  equation  (one-dimensional  Helmholz  equation): 

4r+kz  \EI{z,O))  =  0.  (3) 

\dz  ) 

In  the  air  half-space  (incident  half-space),  the  general  solution  of  this  differential  equation 
can  be  written  as  the  sum  of  the  incident  and  reflected  traveling  waves 

£,(z,ty)  =  £,(  0,ty)[e‘,v  +  f(z,ty)e+,'M],  (4) 

where  f  is  the  reflection  coefficient  of  medium.  In  the  bio-tissue  half-space  (transmission 
half-space),  the  solution  consists  of  only  one  traveling  wave 

Ex(z,(D)  =  T(0,6))-ExMe-ik‘\  (5) 

where  f  is  the  transmission  coefficient. 

If  we  specify  the  electric  field  on  the  plane  z=0  to  be  £,(0,0  and  assume  there  are  no 
reflected  waves  since  we  are  propagating  into  an  infinite  bio-tissue  half-space,  then  using 
the  inverse  Fourier  transform,  Eq.  2,  the  space-time  solution  is 

Ex  (Z,  t)  =  —  P~f(0,  (O)  ■  £,  (0,  (0)ei(to-k:Z)d(0,  z  >  0.  (6) 

2  7tJ-~ 

In  the  following  we  will  make  use  of  the  fact  that  the  incident  pulse  £,(0,0  and  the 
transmitted  pulse  Ex(z,t)  can  be  considered  as  the  input  and  output  of  a  linear  system 
with  transfer  function,  H(z,0J).  Thus,  if ‘the  incident  pulse  £,(0,0)  whose  Fourier 
spectrum  is  £,  (0,  co)  at  z-0  enters  a  linear  medium  then 

Ex(z,0))  =  H(z,  0))Ex  (0,6»,  (7) 

where  £,(z,£0)is  the  Fourier  spectrum  of  the  transmitted  pulse,  and  H{z,cd)  is  the  steady- 
state  transfer  function  of  the  medium.  Eq.  7  describes  the  transmitted  pulse  spectrum 
with  respect  to  depth  z,  the  distance  from  the  planar  interface  separating  the  two  half- 
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spaces.  Eq.  7  is  applicable  when  the  input  and  output  are  both  zero  at  time,  t=0,  which 
will  be  assumed  hereafter. 

Eqs.  1,2  &  7  provide  a  systematic  method  for  finding  the  time-history  of  a  transmitted 
pulse  in  a  linear  dispersive  medium.  Starting  with  the  incident  pulse  Ex(0,t)  obtain  its 

Fourier  transform  from  Eq.  1 


Ex  (0,(0)  =  \2  Ex  (0,  t)emdt .  (8) 

Next,  use  Eq.  7  to  determine  Ex(z,(0)  from  which  Ex(z,t)  is  obtained  using  the  inverse 
Fourier  transform  given  by  Eq.  2 


£*(*.0  =  —-  ]Ex(z,0))e“°'d(0.  (9) 

2  K  J 

Combining  Eqs.  7-9  yields  an  expression  for  the  transmitted  pulse  at  depth  z  and  time  t 

E x(z,t)  =  —  f  H(z,(o)Ex(0,(o)eimd(O,  (10) 

2tt  J-~ 


where:  H(z,cd)  =  f(0,G))e  '*z.  (11) 

The  Fresnel  transmission  coefficient,  f(0,(o),  in  Eq.  1 1  accounts  for  the  effect  on  the 
amplitude  and  phase  of  the  transmitted  pulse  Ex(z,t)  due  to  the  planar  boundary  at  z  = 0 
separating  the  two  half-spaces.  It  is  usually  expressed  in  terms  of  the  indices  of  refraction 
of  the  two  half-spaces.  Assuming  the  incident  half-space  to  be  air  with  index  of  refraction 
nearly  equal  to  unity,  the  Fresnel  transmission  coefficient  at  normal  incidence  is 

f(co)  =  — — - andn(co)  =  vt(to),  (12) 

l+n(G)) 

where  n(cti)  is  the  complex  (frequency  dependent)  index  of  refraction  of  the  half-space 
containing  the  transmitted  pulse  and  l(C0)  is  the  complex  (frequency  dependent)  relative 
permittivity  of  the  medium. 

The  complex  exponential  factor:  exp {-i(kz-(Ot)}  in  Eq.  11  effects  the  propagation  in 
time  and  space  of  the  waveform  Ex(z,t)  by  applying  the  proper  phase  shift  and 
attenuation  to  each  sinusoidal  component  of  Ex(z,(0)  to  produce  a  traveling  wave.  In 
general,  the  wavenumber  it  (to)  =  to  /  v  (to)  and  phase  velocity  are  both  complex  and 
frequency  dependent.  The  real  part  of  the  wavenumber  describes  the  phase  shift  per  unit 
distance  traveled  by  each  sinusoidal  component  while  the  imaginary  part  of  the 
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wavenumber  accounts  for  the  exponential  amplitude  attenuation  per  unit  distance  traveled 
by  each  sinusoid.  The  frequency  response  of  a  dispersive  medium  is  determined  by  the 
frequency  dependence  of  the  complex  permittivity  of  the  medium  and  its  relationship  to 
the  wavenumber  at  frequency  co .  For  bio-tissue,  with  relative  magnetic  permeability  of 
unity,  we  have 


k(a»  =  (a>/c)JWa),  (13) 

where  c  is  the  free-space  velocity  of  light.  The  wavenumber  (or  equivalently  the  relative 
permittivity)  is,  in  turn,  dependent  on  the  medium  in  which  the  wave  propagates. 

To  simplify  the  pulse  propagation  calculations  in  this  report  we  will  make  the  following 
assumptions:  1.)  all  materials  are  nonmagnetic,  2.)  the  dielectric  materials  (bio-tissue)  are 
isotropic  and  homogeneous,  3.)  plane  waves  are  normally  incident  on  a  planar  interface 
separating  air  and  bio-tissue,  and  4.)  the  single  planar  interface  separating  the  two  half¬ 
spaces  is  infinite  in  extent  so  that  diffraction  may  be  ignored.  For  high  frequency 
wideband  pulses  propagating  in  dispersive  tissue,  we  believe  these  assumptions  will  still 
yield  computational  results  accurate  to  first  order  and  have  heuristic  value. 

Debye  Medium 

Liquid  water  is  a  principle  constituent  of  all  living  organisms.  Protoplasm,  the  basic 
material  of  living  cells,  consists  of  a  solution  in  water  of  fats,  carbohydrates,  proteins, 
salts.  Water  acts  as  a  solvent,  transporting,  combining,  and  chemically  breaking  down 
these  substances.  Water  plays  a  key  role  in  the  metabolic  breakdown  of  proteins  and 
carbohydrates.  The  human  body  is  65%  water  by  weight;  some  tissues  such  as  the  brain 
and  lung  are  nearly  80%  water.  It  is  well  known  that  the  dielectric  properties  of  tissue  are 
primarily  determined  by  tissue  water  content  with  both  permittivity  and  conductivity 
increasing  with  increases  in  the  concentration  of  tissue  water  at  RF  and  microwave 
frequencies. 

At  the  microscopic  scale,  the  bonding  of  the  hydrogen  atoms  to  the  highly  electronegative 
oxygen  atom  makes  water  a  polar  molecule.  The  triangular  geometry  of  the  water 
molecule  together  with  its  electrical  polarity  produce  strong  intermolecular  bonds  that 
account  for  water’s  unique  combination  of  physical  properties:  a  high  boiling  point,  a 
solid  phase  that  is  less  dense  than  liquid  pha?e  (ice  floats!),  high  specific  heat,  excellent 
solvent  properties,  and  high  dielectric  permittivity. 

The  triangular  structure  of  the  water  molecule  which  gives  rise  to  a  permanent  dipole 
moment  is  the  source  of  the  polarization  response  of  water.  This  dipole  moment  causes 
the  water  molecules  to  resist  the  random  orientation  due  to  thermal  agitation  and  align 
themselves  with  the  direction  of  an  impressed  microwave  field.  (The  applied  electric 
field,  in  fact,  biases  the  orientation  of  the  molecules  to  only  a  small  extent,  in  comparison 
to  the  thermal  motion  of  the  water  molecules  [7].)  An  a.c.  microwave  field  will  produce  a 
sinusoidal  torque  that  will  cause  individual  water  molecules  or  molecular  groups  to  rotate 
with  an  angular  velocity  proportional  to  the  torque.  The  angular  momentum  of  the  water 
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molecules  results  in  friction  with  neighboring  molecules  and  converts  thereby  to  linear 
momentum,  which  by  definition  is  heat  in  liquids  and  gases.  This  response  of  the  water 
molecules  is  called  relaxation  polarization  and  is  dependent  on  the  frequency  of  the 
electric  field.  P.  Debye  [6]  was  the  first  to  successfully  model  the  relaxation  polarization 
response  of  water  as  due  to  the  rotational  motion  of  the  water  molecule  in  a  damping 
frictional  medium.  See,  for  example,  references  [7,8]  for  detailed  accounts  of  modem 
theories  of  the  macroscopic  dielectric  response  of  water  given  in  terms  of  molecular 
quantities. 

The  Debye  equation  has  been  used  to  model  the  frequency-dependent  polarization 
response  of  water  and  ice  to  an  impressed  microwave  field  over  the  frequency  range  0.3 
to  300  GHz.  The  complex  relative  permittivity  of  water  modeled  by  the  Debye  equation 

is 


m  =  em+  .  <14> 

1  + 1  COT  COBq 

where  z  is  the  dielectric  relaxation  time,  and  0)  is  the  angular  frequency  of  the  impressed 
time-varying  field,  and  £0  =  10-9  136k  (F/  m)is  the  permittivity  of  free-space.  The 

relaxation  constant,  r  =  8.1  x  10“12  s ,  measures  the  decay  time  of  the  macroscopic 
polarization  of  water  when  the  external  field  is  removed.  It  corresponds  to  the  water 
absorption  peak  located  near  20  GHz  from  the  relation  on=l.  The  constants 
e  =  55  and  es  =  78.2  are  the  limiting  values  of  the  effective  dielectric  constants  of  water 
at  the  high  and  low  frequencies,  respectively.  The  ionic  conductivity  of  pure  water  at  25 
°C  is  represented  by  the  constant  a  =  10-4  Siemans  /  m .  These  numerical  values  are 
temperature  dependent  and  were  obtained  by  empirical  curve  fitting  methods  applied  to 
dielectric  measurements  [2].  (Note:  the  departure  of  £„,  from  unity  is  due  to  the  fact  that 
only  a  limited  amount  of  data  spanning  a  finite  range  of  frequencies  were  in  used  to 
estimate  the  parameters  in  the  Debye  equation.)  These  parameter  values  will  be  used  in 
the  following  computations  for  water  based  tissue  described  as  a  single  relaxation  time- 
constant  Debye  medium. 

Lorentz  Medium 

Another  model  for  tissue  absorption  and- dispersion  in  the  range  3  x  1014  to  3  x  1015  Hz 
(i.e.  in  the  optical  range  from  the  u.v.  (0.1  |im)  to  the  near  i.r.  (1  fim) )  comes  from  the 
mechanical  model  of  electron  resonance  developed  by  Lorentz  and  Drude.  This  model  of 
electronic  resonance  is  not  exact,  since  it  ignores  the  fact  that  the  electron  s  motion 
follows  the  laws  of  quantum  mechanics  rather  than  newtonion  mechanics. 
Notwithstanding  the  known  limitations  inherent  in  the  mechanical  analogy,  this  model 
gives  good  results  [5].  This  dielectric  model  describes  the  resonance  polarization  in  a 
homogeneous,  isotropic  dielectric  characterized  by  a  single  resonance  frequency  co0  for 

which  the  complex  relative  permittivity  is  given  by 
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(15) 


.  1 
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Here,  b2  is  the  square  of  the  plasma  frequency  of  the  medium,  6  is  a  damping  constant  to 
account  for  the  energy  lost  by  the  harmonically2  bound  electron  at  the  resonant  frequency 
The  damping  that  occurs  is  not  due  to  the  electron  moving  in  a  viscous  fluid,  but  is 
due  to  the  energy  lost  by  the  electron  by  radiation  as  the  result  of  forced  oscillations. 

The  medium  parameters:  co0  =  4  x  10^  Is,  =  20 x  10 ^  ,  and  S  =  0.28 x  10^  /  s 

describe  a  strongly  dispersive  and  absorptive  medium  with  an  absorption  band  in  the 
ultraviolet  region  of  the  electromagnetic  spectrum.  These  medium  parameters  were 
suggested  in  [3]  and  are  the  same  numerical  values  used  by  Brillouin  in  [9,10].  They  will 
be  used  in  the  following  computations  for  the  single  resonance  Lorentz  medium. 


Pulse  Types 

It  remains  to  define  the  two  types  of  transient  waveforms  to  be  discussed.  The  specific 
waveforms  addressed  in  this  report  are  (1)  the  finite  train  of  pulsed  sinewaves  and  (2)  the 
single  Gaussian  modulated  sinewave  pulse.  These  two  waveforms  are  prototypical,  if 
somewhat  idealized,  of  those  generated  by  real  pulsed  radio  frequency  and  laser  emitters. 
Perhaps  more  importantly  for  this  discussion,  their  Fourier  spectra  can  be  found 
analytically. 

Sinewave  Pulse 

With  the  aid  of  the  Heaviside  unit  step  function,  H(t), 


H(t-a)  = 


t<a 
t>a  ' 


(16) 


we  can  write  a  concise  analytical  expression  for  any  finite  pulse  train. 

We  consider  first  a  single  pulse  of  sinewave  carrier  ©..amplitude  A,  duration  i  at  z=0  as 
illustrated  in  Figure  1.  This  can  be  written  as 


f(t)  =  A-sin((Oct) 


V  2 


J 


(17) 


Its  Fourier  transform  is  pure  imaginary,  is  derived  in  Appendix-A,  (see  Eq.  A-9)  and  can 
be  written  as  follows: 

F(a>)  =  -  (sinc(or)  -  sinc(/l)) ,  (18) 

where  i  =  V-l;  sinc(x)  =  sin x / x;  a  =  ( (Oc  -q))t/ 2;  fi  =  ((Oc  +  (0)T/2. 


2  the  restoring  force  is  linearly  proportional  to  displacement 
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The  spectrum  of  this  sinewave  pulse  is  similar  to  the  spectrum  of  its  rectangular  envelope 
but  translated  in  frequency  by  an  amount  equal  to  the  carrier,  mc.  The  width  of  the  main 
frequency  lobe  increases  as  the  length  of  the  rectangular  modulation  waveform,  z, 
decreases  and  is  independent  of  the  carrier  frequency.  The  first  sine  term  accounts  for 
most  of  the  energy  located  near  the  carrier  frequency,  while  the  second  sine  term  becomes 
significant  only  at  frequencies  far  removed  from  the  carrier.  It  follows  also,  that  the  value 
of  |F(<u)|  at  to  =  toc  is  |F(<a)|  «  At/2.  In  Figure  2  the  magnitude  of  the  spectral  density 
function  of  a  single  sinewave  pulse  of  duration  T  is  plotted  over  the  positive  frequency 
axis.  We  see  that  the  mainlobe  peak  is  centered  at  the  carrier  frequency  cot.  and  has 
spectral  width  of  Aco  =  2  /  z .  It  is  interesting  to  note  that  when  this  pulse  contains  an 
integer  number,  N,  of  carrier  cycles,  i.e.  z  =  2nN  /  o)c ,  then  the  number  of  zero-crossings 
in  the  spectral  amplitude  function  between  zero  and  <ac  is  also  equal  to  N.  This  is  true 
because  sinc(a)  =  sinc[^V(l -tolcoc)]  =  0;  when  co / 0)c  =  1  /  N,  2 /  N . . ,(N - 1)  /  N. 

We  consider  next  a  finite  number,  M,  of  sinewave  pulses,  where  M>1.  Without  loss  of 
generality,  we  take  the  number  of  pulses  in  the  finite  pulse  train  to  be  odd:  ( M=2N-1 ), 
where  the  middle  pulse  is  centered  at  the  origin,  the  pulse  width  is  Z  ,  and  the  pulse 
period  is  t0 .  The  finite  pulse  train  consisting  of  an  odd  number  of  pulses  can  be  written  as 

follows 


N.(  r 

H 

„=i  V  . 


f(t)  =  Asin(fi)ct)X 


+Asin(to(.t) 


+Asin(<0t.O^  H 

n=l  V 


t  +  ntn  — 
0  2 


-H 


t+ntn  H — 


( 

T 

-H 

r  *ti 

H 

t  +  — 

\ 

L  2j 

L  2j  J 

(19) 


X 

-H 

tT 

t-ntn  -l — 

t-nt0—— 

2j 

2j  J 

(Note:  Eq.  19  is  equivalent  to  Eq.  A-l  in  Appendix-A).  This  finite  duration  wavetrain  is 
generated  by  gating  or  “turning  on”  each  pulse  at  a  positive  zero-crossing  for  an  integral 
number  (nO  of  carrier  cycles,  followed  by  an  integral  number  (n2>  of  “off ’-cycles, 
thereby  creating  a  coherent  pulse  train  with  pulse  repetition  period  equal  to  (m+n2)  carrier 
cycles.  A  wideband  pulse  will  be  created  when  only  a  few  carrier  cycles  (m)  occur  during 
the  pulse  “on”  time.  The  envelope  of  the  pulse  train  is  amplitude  modulated  using  a 
rectangular  gating  function  requiring  maximum  bandwidth  in  the  Fourier  spectrum.  The 
Fourier  transform  of  a  unit  amplitude  pulse  train  is  derived  in  Appendix-A  and  is  shown 
to  be 
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2i  sin[(fflr -nm0)  Afa  /a  Jsin[(ft)  +na0)772] 

L/rU)J  -  ^  (fflc  -nG)0)(m  +  n(W0) 

2i  «y  sin[(6)c  +  nw^Nn  /  a)c  ]  sin[(ft) + mo0  )T  /  2] 

7q  n=-«  (U)c  -t~ wO)q)(U) -t* H(Oq) 

(20) 

_  2/  «y°  sin[(6)c  -  »6)0)T  /  2]  sin[(6) + nft)0)r /  2] 

_f0„=^  (©c-nft)0)(6)  +  n<O0) 

2i n  =  sin[(ft)c  +n6)0)T/2]sin[(<Q+»6?0)r/2] 

r0  „=—  (coc+nQ}0)(a)+nffl0) 

For  a  finite  duration  wavetrain  consisting  of  M  pulses,  the  spectrum  will  consist  of  a 
series  of  lobes  of  width  1  /  Mt0  modulated  by  an  envelope  which  corresponds  to  the 
spectrum  of  one  pulse.  The  separation  between  lobes  is  1  /t0  =  PRF . 

Gaussian  Pulse 

A  bell-shaped  pulse  can  be  obtained  by  impressing  a  Gaussian  amplitude  modulation  on  a 
sinewave  carrier.  At  /  =  0  this  waveform  may  be  expressed  as  the  real  part  of 

f{f)  =  Ae-'2/(2f)V'<8"',  (21) 

where  V2 Z  =  standard  deviation.  Here,  A  is  the  peak  amplitude, &)0  is  the  angular 
frequency  of  the  sinewave  carrier  and  %  is  the  pulse  width  parameter.  Since  a  time  span 
of  four  standard  deviations  centered  at  t  =  0  contains  approximately  95%  of  the  area 
under  the  Gaussian  curve,  a  wideband  pulse  will  therefore  exist  when  only  a  few  periods 
of  the  sinewave  carrier  occur  during  this  time  span.  The  Fourier  transform  of  this  pulse  is 
derived  in  Appendix-B. 


Numerical  Procedures 

Positive  Frequencies  Only 

To  find  the  time-history  of  the  propagated  waveform  for  any  time  t  >  0,  it  is  necessary  to 
evaluate  Eq.  10  by  numerical  quadrature  since  a  general  analytic  solution  is  unavailable. 
However,  we  will  show  that  it  is  only  necessary  to  evaluate  the  inverse  Fourier  transform 
over  positive  frequencies  by  using  a  well-known  property  of  real-valued  signals.  For  real 
signals,  the  complex  conjugate  of  Eq.  1  is  given  by 


F*m=]nt)ei<0,dt. 


(22) 


From  Eq.  22  follows  the  Hermitian  symmetry  of  the  Fourier  transform  for  real-valued 
signals  so  that 
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V. 


F(-o))  =  F*(o)).  (23) 

This  property  of  the  Fourier  transform  of  real-valued  signals  implies  that 


|F(-ty)|  =  |F*(fi))|=|F(ty)|.  (24) 

Whatever  shape  the  modulus  of  the  spectral  density  function  has  for  positive  frequencies, 
it  is  the  same  for  negative  frequencies,  i.e.  only  one-half  of  the  real  signal  spectrum  is 
needed  to  uniquely  specify  the  signal.  We  divide  the  integral  in  Eq.  2  into  two  parts  as 

/(f)  =  —  J  F(co)eim'dco+—]  F(0))e“°'d0).  (25) 

2ft  2k  0 

By  use  of  Eq.  23  we  can  write 

f(t)  =  —][F(co)em,+F\o))e-im]dco.  (26) 

2k  0 

We  note  that  the  quantity  in  brackets  is  real.  Therefore,  we  can  rewrite  Eq.  26  as 


2k  •  f(t)  =  Rejj  2F(0))eiw,do) 

This  shows  that  it  is  only  necessary  to  perform  the  numerical  quadrature  over  positive 
frequencies  to  recover  the  propagated  real  signal. 

Error  Analysis 

In  order  to  determine/ft)  in  Eq.  2  or  Ex(z,t)  in  Eq.  10,  integration  must  be  performed 
along  the  <3J-axis.  To  perform  the  integration  we  choose  to  uniformly  sample  the 
integrand  at  a  finite  resolution  Act)  over  a  finite  frequency  range  Q .  We  consider  the 
errors  resulting  from  each  of  these  finite  approximations  and  consider  the  errors  due  to 
sampling  first. 

Aliasing  Error 

Implementing  the  numerical  quadrature  of  the  inverse  Fourier  transform  by  means  of  the 
rectangular  rule,  with  fixed  step  size,  Aco,  we  can  approximate  Eq.  2  by  writing  it  as  an 
infinite  sum  over  the  positive  integers  which  corresponds  to  folding  the  integral  over 
positive  frequencies  as  shown  in  Eq.  27 

fs (t)  = — Y  F(nAo})einAat .  (28) 

ft  n= o 


10 


We  recognize  Eq.  28  as  a  complex  Fourier  series  representation  of  fs(t).  To  make  this 
explicit,  we  invoke  Eq.  A- 12  from  Appendix-A  and  let  Ct)0  =  Act)  and 
F(nAco)  =  cnt0  =  cn2n/  A  0)  to  obtain 


fs(t)  =  2YjCne“ 


mAcor 


n- 0 


(29) 


By  uniformly  sampling  the  continuous  Fourier  spectral  density  in  order  to  numerically 
approximate  the  Fourier  integral,  we  arrive  at  a  Fourier  series  representation  of  the 
function  fft)  with  period  t0  =  2nl  Act) . 

We  emphasize  that  the  recovered  time  function  fftfis  periodic  with  repetition  period 
given  by  t0  =  2n  /  Ace.  This  can  also  be  demonstrated  by  replacing  t  in  Eq.  28  with 
t  +  m-2nl  Act)  where  m  is  any  integer.  We  obtain 

f,(l+m-2lc/Ao»  = — jr 

ft  n= 0 

= — V  F{nAo))einAea+inm{2!t)  (30) 

^  n=0 

=— Y,F{nA(0)einAc* 

ft  n= 0 

Thus  we  see  that  fs(t  +  m-2n/  Aco)  =  fs(t),  i.e.  a  function  sampled  in  the  frequency 
domain  will  be  repetitive  in  the  time  domain  with  period,  t0  =  2n  /  Aco.  We  must  choose 
Act)  small  enough  so  that  tQ  is  larger  than  a  time  interval  which  contains  most  of  the 
energy  in  the  signal,  otherwise  an  unacceptable  amount  of  aliasing  error  will  occur  in  the 
reconstructed  time  domain  signal. 

Removing  the  Periodicity 

The  disadvantage  of  the  Fourier-series  representation  of  Eq.  28  lies  in  the  fact  that  fs(t)  is 
a  periodic  function,  so  it  cannot  be  used  io.approximate  the  aperiodic  time  domain 
waveform  directly.  However,  it  is  possible  to  retain  only  that  part  of  fft)  inside  the 
interval  0  <  t  <  t0  to  recover  an  aperiodic  waveform  [11].  To  accomplish  this,  we  first 
multiply  1ft)  by  the  unit  step  function  to  remove  that  part  of  fft)  in  the  range  t  <  0 .  To 
eliminate  the  part  of  fft)  for  t  >  t0,  we  can  use  the  double-impulse  weighting  function 
w(t)  given  by 

w(t)  =  S(t)-d(t-t0).  (31) 
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Convolving  Eq.  31  with  /,(/)  is  equivalent  to  adding  to  fs(t)  its  own  negative  delayed  by 
one  period.  So  performing  these  two  operations  in  the  complex  frequency  domain,  by 
means  of  the  Laplace  transform  we  have 


n= 0  * 


2F(nAcd)  /  f0 


s—inAo) 


(32) 


which  leads  to  the  desired  aperiodic  time-domain  representation 


71=0  *0 


where //(0  is  the  unit  step  function  defined  by  Eq.  16.  We  observe  that  the  first  infinite 
sum  of  Eq.  32  is  the  Laplace  transform  of  the  semi-infinite  periodic  function  fs(t)  starting 
at  t  =  0  and  the  second  infinite  sum  of  Eq.  32  is  the  Laplace  transform  of  an  identical 
semi-infinite  periodic  function,  fx(t  —  t0),  starting  at  t  =  t0  that  exactly  cancels  the  first 
function  from  t  =  tQtot  =  °°. 


Also,  F(z,nAco)  =  Ex{z,co )  =  H{z,co)  ■  £,(0,<a )  from  Eq.  3,  where  H(z,co)  is  the 
medium’s  transfer  function  (not  to  be  confused  with  H(t),  the  Heaviside  step  function) 
and  Ex(z,co)  is  the  Fourier  transform  of  the  propagating  pulse  at  depth,  z.  If  the  time 
domain  function  £,(0,0  has  duration  7J  and  the  “filter”  impulse  response  H(z,t),  has  a 
duration  T2 ,  then  the  convolution  result: 

Ex(x,t)  =  H(z,t)®Ex(0,t),  (34) 


has  a  duration  equal  to  Tx  +T2.  If  the  “filtering”  is  done  by  multiplication  in  the  frequency 
domain, 

F(z,nAco)  =  Ex{z,co)  —  H(z,co)  •  E x(0,a>) ,  (35) 


and  Aa  is  chosen  as  In  / 10  with  t0<(Tx+T2),  then  overlapping  in  the  time  domain  will 

2  71 

occur  and  Ex(z,t)  will  be  distorted.  If  7J  >  T2,  then  Ao)  =  —  is  a  safe  choice  and  the 

V.  ^*1 

overlapping  effect  (aliasing)  in  the  time  domain  cannot  occur. 

Truncation  Error 

Next,  we  consider  the  error  resulting  from  the  evaluation  of  the  integral  over  a  finite 
frequency  range,  Cl .  By  truncating  the  integrand  in  Eq.  2  we  obtain  a, finite  inverse 
Fourier  transform 


-  1  ft/2 

/(f)  \F(co)eim,dco. 
2ft  -an 


(36) 
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Because  the  spectral  densities  of  the  pulsed  waveforms  under  discussion  have  infinite 
tails,  it  is  necessary  to  examine  the  relationship  between  the  error  in  the  computed  inverse 
transform  and  the  finite  frequency  interval  Q .  By  knowing  this  relationship  it  is  possible 
to  select  the  proper  limits  on  the  inverse  transform  to  achieve  a  desired  accuracy. 

The  finite  inverse  Fourier  transform  can  be  written  as  the  infinite  transform  of  a  product 
of  functions  when  one  of  those  functions  is  the  rectangular  window  function, 

B{(o)  between  Q.  /  2  and  -Q./2] 


where 


f{t)  =  -^\l'F{co)B{(D)eicad(0, 


(37) 


B(a»  = 


-£2/2<co  <Q/2 
elsewhere 


(38) 


The  inverse  of  the  rectangular  window  function,  B(co),  is  obtained  as  follows 


2k 


an 


-an 


1 

\e'nn' 

(2/);rr 

L 

1  •  1 

=  — sin 

— 

Kt  \ 

<2) 

Q  . 

(Clr 

=  — sine 
2k 


(39) 


Multiplication  in  the  frequency  domain  is  equivalent  to  convolution  in  the  time  domain. 
Therefore,  we  have  for  a  pulse  centered  at  the  origin  with  duration  equal  to  r0 


+t0/2 

/( 0=  jf(t-T)b(T)dr.  (40) 

-T0/2 

For  example,  if  f(t)  were  a  single  sinewavb  pulse  of  unity  amplitude  and  carrier  o)c,  we 
would  have 

f(t)= —  Jsin[toc(r— T)]-[H(r-T  +  T0/2)  — if(f— T-r0/2)]-sinc|^— Jrfr.  (41) 

-T0/2 

Eq.  40  shows  clearly  how  the  finite  inverse  Fourier  transform  affects  the  computed  time 
domain  pulse.  The  distorted  waveform  is  obtained  by  convolving  the  true  value  of  the 
pulse  waveform  with  the  sine  function  for  each  time  point.  The  difference,  f(t)-f(t),  is 
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the  error  due  to  the  finite  integration  limits.  The  error  obtained  in  this  way  is  independent 
of  the  aliasing  error  due  to  sampling  as  discussed  previously.  Figures  16-18  show 
examples  of  this  effect  on  a  1  radian/sec,  four  cycle,  unit  amplitude  sinewave  pulse.  The 
distortion  of  the  pulse  caused  by  evaluating  its  inverse  Fourier  transform  over 
progressively  smaller  bandwidths  of  Q  at  10,  5  and  2.5  radians/sec  is  demonstrated. 

Results  and  Discussion 
Sinewave  Pulse  in  a  Debye  Medium 

Our  first  set  of  results  models  the  propagation  of  a  single  sinewave  pulse  with  carrier 
frequency  1  GHz  normally  incident  onto  a  water  half-space  (Debye  medium).  This 
incident  pulse,  shown  in  Figure  1,  is  a  single  burst  of  10  carrier  cycles.  Figure  3  shows  the 
transmitted  pulse  plotted  against  time  at  a  depth  of  0.75  meters  in  water.  We  note  the 
attenuation  with  depth  and  the  presence  of  the  Brillouin  precursors  at  the  leading  and 
trailing  pulse  edges.  This  result  duplicates  that  obtained  by  previous  investigators  [1]  who 
used  a  Fourier  series  calculation.  The  FORTRAN  code  listing  for  this  case  is  included  in 
Appendix-C. 

Figure  4  illustrates  the  propagation  of  2  1/2  pulses  of  a  5  pulse  wavetrain  in  a  water  half¬ 
space  showing  how  rapidly  the  carrier  is  attenuated.  At  1  meter  depth  all  that  remains  of 
the  pulse  is  the  Brillouin  precursor  pair  associated  with  the  leading  and  trailing  edges  of 
the  pulse.  Figure  5  shows  the  amplitude  spectrum  plotted  over  the  positive  frequencies  for 
a  finite  duration  pulse  train  consisting  of  five  such  pulses.  Compared  to  the  amplitude 
spectrum  of  a  single  pulse  shown  in  Figure  2,  Figure  5  reveals  a  more  complex  lobe 
structure  due  to  the  infinite  sum  of  sine  functions  needed  to  describe  the  spectrum  of  the 
five  pulse  train.  Eq.  A- 19  and  subsequent  discussion  in  Appendix-A  explain  the  origin  of 
this  more  complex  structure. 


Sinewave  Pulse  in  a  Lorentz  Medium 

Our  second  set  of  results  investigates  the  propagation  of  a  finite  sinewave  pulse  train  (five 
pulses)  normally  incident  onto  a  Lorentz  half-space  with  a  single  resonance.  As 
mentioned  previously,  the  Lorentz  model  is  appropriate  for  describing  the  absorption  and 
transmission  of  infrared  to  optical  signals  in  certain  materials.  Figure  6  shows  two 
“snapshots”  of  a  five-pulse  burst  taken  at  6*0  and  100  femtoseconds  after  incidence.  The 
carrier  frequency  of  this  pulse  train  lies  near  the  upper  end  of  the  absorption  band.  The 
band-stop  filtering  action  of  the  Lorentz  medium  blocks  the  pulse  energy  that  falls  within 
its  absorption  band.  Only  the  energy  outside  the  absorption  band  is  transmitted.  The 
spectral  residue  falling  below  the  absorption  band  form  a  Brillouin  precursor  pair  for  each 
incident  pulse.  The  spectral  residue  above  the  absorption  band  form  a  Sommerfeld 
precursor  pair  for  each  incident  pulse.  Since  these  two  precursor  types  are  separated  in 
frequency  by  the  width  of  the  absorption  band,  one  would  expect  them  to  move  at 
different  speeds  in  a  linear  dispersive  medium  such  as  this  one.  Comparing  Figures  6  (a) 
&  (b)  indeed  shows  this  to  be  true.  We  see  the  distance  between  the  Sommerfeld  and 
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Brillouin  precursor  pairs  increasing  with  time  which  indicates  a  difference  in  propagation 
speeds.  A  series  of  plots  similar  to  Figure  7  was  generated  to  measure  the  difference  in 
travel  times  for  various  distances  up  to  200  microns.  A  linear  regression  of  displacement 
versus  time  for  the  Brillouin  precursor  pair  in  Figure  8  shows  that  this  pair  propagates  at 
two-thirds  the  speed  of  light.  Similarly,  a  linear  regression  was  done  for  the  Sommerfeld 
precursor  pair  of  displacement  versus  time.  Figure  9  indicates  that  the  Sommerfeld 
precursor  pair  propagates  at  nearly  vacuum  light  speed  in  this  medium.  Using  the  results 
for  the  Brillouin  precursor  speed  and  Snell’s  law  yield  a  refraction  angle  for  the  Brillouin 
precursor  of  about  41°  for  light  obliquely  incident  at  80°  on  a  Lorentz  half-space  which 
agrees  with  the  result  reported  in  [4]. 


Gaussian  Pulse  in  a  Lorentz  Medium 

Our  third  set  of  results  addresses  the  propagation  of  a  wideband  Gaussian  transient  into 
the  same  single  resonance  Lorentz  medium  as  discussed  above.  This  pulse  is  assumed  to 
be  normally  incident  on  the  Lorentz  half-space,  with  pulse  width  parameter,  t,  equal  to 
.05  femtoseconds,  and  carrier  frequency,  (Oc  =  5.75E  +  16  rad  /  s .  The  goal  was  to  match 
the  results  shown  in  Figure  8(a)  of  reference  [3].  Our  computation  is  based  on  the  Fourier 
transform  of  a  cosine  carrier  modulated  with  a  Gaussian  envelope  as  given  in  Appendix- 
B. 

Figure  10  shows  the  spectrum  of  this  pulse.  We  show  in  Figure  1 1  the  time-history  of  this 
pulse  at  a  depth  of  one  micron  in  the  Lorentz  half-space.  We  note  the  close  resemblance 
to  Figure  8(a)  in  reference  [3].  This  result  validates  our  Fourier  integral  numerical  scheme 
by  favorable  comparison  to  another  independent  computational  method  based  on  the 
numerical  inversion  of  the  Laplace  transform.  The  FORTRAN  code  for  this  case  is 
included  in  Appendix-D. 

Impulse  Response  of  a  Dielectric  Medium 

By  setting  Ex(0,0))  equal  to  unity  in  Eq.  10,  we  obtain  the  impulse  response  of  the  linear 
medium 

h(ZJ)  =  -L  r  H(z,Q))eic0,dCQ.  (42) 

2  K 

^  #  . 

The  above  integral  can  be  computed  by  numerical  quadrature  using  the  same  techniques 
previously  discussed. 

The  time-domain  convolution  integral,  provides  an  alternative  means  of  obtaining  the 
propagated  time-domain  pulse  at  depth  z  of  the  incident  waveform £^(0,0  when  the 
impulse  response  of  the  medium  is  available  as  follows 

Ex(z,t)  =  \‘QEx(0,t)h(z,t-T)dT.  (43) 
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In  the  evaluation  of  the  above  convolution  integral  by  digital  computer,  it  is  necessary  to 
replace  the  integral  with  a  finite  sum.  We  first  replace  the  infinitesimal  dz  with  a  time 
interval  of  finite  duration  T,  and  let  nT  =  Z .  Then  we  can  rewrite  Eq.  43  as  the 
summation 


Ex(z,t)  =  TfjEx(0,nT)h(z,t-nT)  for  KT <t  <{K  +  \)T .  (44) 

n=0 

Figure  12  shows  the  impulse  response  of  pure  water  (Grant’s  model)  at  a  depth  of  0.75 
meters  where  h(z,t)  is  the  water  response  measured  at  depth  z  for  an  impulse  incident  on 
the  z  =  0  plane.  In  Figure  13  we  computed  the  time-history  of  a  10  ns  pulse  of  1  GHz  at 
this  depth  by  means  of  Eq.  44.  We  note  the  similarity  to  Figure  3  which  was  computed 
using  the  Fourier  transform.  The  FORTRAN  code  for  this  case  is  included  in  Appendix- 
E. 

Figure  14  shows  the  impulse  response  of  the  single  absorption  Lorentz  medium  at  1 
micron  depth.  Once  again  using  Eq.  44,  we  compute  the  time-history  of  this  pulse  in  the 
Lorentz  medium  and  see  close  agreement  with  Figure  7(a).  The  FORTRAN  code  for  this 
case  is  given  in  Appendix-F. 


Conclusion 

We  have  shown  how  to  compute  the  waveform  of  a  transmitted  pulse  in  a  linear 
dispersive  half-space  by  use  of  the  pulse’s  Fourier  transform.  The  approach  discussed  in 
this  paper  is  applicable  only  to  waveforms  whose  Fourier  transform  can  be  found 
analytically.3  In  particular,  we  have  demonstrated  this  method  for  a  finite  train  of 
sinewave  pulses  and  a  single  Gaussian  modulated  sinewave  transient.  These  two  idealized 
waveforms  reasonably  approximate  those  often  encountered  in  dosimetry  computations 
involving  “thick”  homogenous  layers  (half-space)  of  bio-tissue  and  therefore  should  be 
useful  in  estimating  electric  field  strengths  in  bio-dielectrics  that  can  be  accurately 
described  by  either  the  Debye  or  Lorentz  models. 

An  unavoidable  side-effect  of  uniformly. sampling  the  spectral  density  function  in  order  to 
numerically  compute  the  inverse  Fourier  transform  is  that  the  resulting  time-domain 
waveform  is  periodic  with  period  equal  to  the  reciprocal  of  the  frequency  sampling 
interval.  We  have  shown  how  to  analytically  remove  this  periodicity  and  thereby  recover 
the  transient  waveform. 


3  The  Fourier  transform  method  itself  is  not  limited  to  analytically  describable  waveforms.  It  is  possible  to 
start  with  an  arbitrary  time  waveform  sampled  at  discrete  time  points  and  compute  its  Fourier  transform 
using  the  DFT;  multiply  by  the  discretized  medium  transfer  function  and  invert  the  result  via  the  IDFT  to 
obtain  the  propagated  time-domain  wave  in  the  medium.  The  accuracy  of  this  method  is  limited  by  the 
sampling  resolution  of  the  initial  waveform.  The  desire  to  avoid  this  limitation  was  the  motivation  for 
choosing  only  analytically  describable  signals  in  this  report. 
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Recommendations 

Additional  work  might  profitably  be  spent  on  refining  these  analytically  convenient 
models  of  pulsed  waveforms  to  make  them  more  closely  resemble  real-world  pulses.  For 
example,  the  sinewave  pulse  should  be  modified  to  have  a  nonzero  rise-time,  as  this  is 
true  for  all  band-limited  pulses  of  physical  origin.  On  the  other  hand,  the  Gaussian  pulse 
should  be  modified  to  have  truncated  tails  which  is  also  true  for  all  physical  pulses  of 
finite  duration. 

To  increase  the  utility  of  these  numerical  pulse  propagation  models,  it  would  be  desirable 
to  speed-up  the  computations  by  harnessing  the  speed  of  the  Fast  Fourier  Transform 
(FFT)  algorithm  to  evaluate  the  inverse  Fourier  transform.  That  is,  the  summation 
inherent  in  the  evaluation  of  the  Fourier  integral  by  means  of  numerical  quadrature  would 
be  performed  by  the  inverse  FFT  algorithm. 
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Appendix  A  Pulse  train  spectrum 

All  time  domain  pulsed  signals  can  be  described  in  terms  of  four  time  scales:  wave 
period,  pulsewidth,  pulse  repetition  interval  (PRI),  and  dwell  time.  Likewise,  in  the 
frequency-domain,  pulsed  signals  can  be  described  in  terms  of  four  corresponding 
frequency  scales:  carrier  frequency,  pulse  bandwidth,  pulse  repetition  frequency  (PRF), 
and  spectral  line  bandwidth  [Al].  In  Fourier  analysis  an  inverse  relationship  exists 
between  the  time  and  frequency  scales:  the  larger  the  wave  period,  the  smaller  the  wave 
frequency. 

This  derivation  of  the  Fourier  transform  of  a  finite  number  of  sinewave  pulses  will 
illustrate  the  four  frequency  scales  mentioned  above. 

For  the  finite  regularly  repeating  signal  shown  in  Fig.  Al,  fr{t)  is  the  product  of  the 
infinite  duration  periodic  function  /(/),  and  the  gate  function  g(t) 

fr(t)  =  f(0- g(t)  (A'l) 

where 


and 


/  (t)  =  sin(6Jcr), 
=  0, 

=  sin(<ycr), 


-?/2<r<7/2 
r/2<t<tQ-r/2 
-t0  +  r/2<t  <-t/2 
t0-T/2<t<t0  +  T/2 
-t0-T/2<t<-t0  +  r/2 


=  0, 


[  t0  +  r/2<t<2t0-r!2 
\-2tn+r/2<t<-t()-T/2 


etc.. 


g(t)  = 


(1,  -TI2<t<T/2 
1 0,  elsewhere 


(A-2) 


(A-3) 


Let  an  integer  number,  N, 
t  =  2nN  I  (0 c. 


of  carrier  cycles  occur  during  the  pulse  “on-time”,  t ,  such  that 


We  begin  by  deriving  the  Fourier  transform  F0(<b)  of  a  single  sinewave  pulse,  of  finite 
duration,  t ,  unit  amplitude,  and  symmetric  with  respect  to  the  origin,  given  by  the  first 
line  of  Eq.  A-2,  which  we  call  f0(t) .  We  use  the  Fourier  transform  pair  definition 


eo  1 

F{(Q )  =  J  f(t)e~iaIdt  and  fit)  =  —  J  Fia)eimd(0 


(A-4) 
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where  i  =  V-l. 

Since  f0(t)  is  a  finite  duration  sinewave  and  the  sinewave  is  an  odd  function  we  see  that 


F[f0(t)]  =  F0(co)  =  -i  j  f0(t)  sin  (cot)dt 


(A-5) 


and 


1/2 

F0(co )  =  -i  J sin(tocr)  sin (cot)dt 


(A-6) 


-t/2 


Note  that  since  f0(t) is  real  and  odd,  F0(a)will  be  imaginary  and  odd.  Making  use  of  a 
trigonometric  identity  we  get 


F0(co)  =  -i  J 1  /  2{cos[(tyc  -  (o)t ]  -  cos[(t»c  + 

-t/2 

Performing  the  integration,  we  have 
F0(o))  =  -i< 

_  .  ,fsin[(ty/.-<y)T/2]  sin[(<yf  +©)r/2]] 

o(fy)  =  'li  (®e  +  «)  J 


sin  [{coc-(o)t] 

T/2 

sin[(tyt.  +  o))t] 

r/2  ' 

V. 

2(6) c- CO) 

-t/2 

2  (6)c+6)) 

-Til  j 

(A-7) 


(A-8) 


(A-9) 


.fsin[(ty  -co)Njtl a)c]  sin[(6>c  +  6?)N^~/6Jc] 
Fo(")_  l<  (c oc-a ))  ( coc  +  a >) 


(A-10) 


Now  it  can  be  shown  [A2],  that  the  complex  coefficients,  cn ,  of  the  Fourier  series 
expansion  of  the  periodic  pulse  train  of  sinusoidal  pulses,  f(t),  with  repetition  period,  f0, 
are  equal  to  the  values  of  the  Fourier  transform  F0(<a)  of  the  single  sinewave  pulse,/,  (0, 
evaluated  atro  =  nco0  =  nln  / 10  and  multiplied  by  1  / 10,  where/,  (r)is  defined  by 


J/o (0.  l'N'o/2 

[0,  \t\>t0/2 


(A-ll) 


Therefore,  the  Fourier  series  coefficients  of  the  infinite  duration  periodic  pulsed  sinewave 
train,  f(t),  defined  by  Eq.  A-2,  are  given  by 
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(A-12) 


—  F0(nco0) 


where 


i  fsin[(ft),. -nft)0)Afrr/ft)c]  sin[(fl)c +tttP0)Mr/fl)J 
t0\  ( (Oc-nCO0 )  (£Oc+nf»0) 


/(0  =  flXcncfa‘B»';  6)0=2x/ t0  (A-13) 

It  can  also  be  shown  [A2]  that  the  Fourier  transform  of  any  infinite  train  of  pulses  is  given 
by 


rt=+«o 

F{co)  =  2n'^cn  S(ot)-no)0 ) 


(A-14) 


where  cn  are  the  complex  coefficients  of  the  Fourier  series  expansion,  and  6  is  the  unit 
impulse  function. 

That  is,  the  Fourier  transform  of  a  periodic  function  consists  of  a  sequence  of  equidistant 
impulses  located  at  the  harmonic  frequencies  of  the  function  with  amplitudes  equal  to  the 
complex  coefficients  of  its  Fourier  series  expansion.  To  obtain  the  Fourier  transform  of 
th  &  finite  pulse  train  of  duration  T,  we  need  the  Fourier  transform  of  the  gating 
function,  g(t) ,  given  by  Eq.  A-3.  Its  Fourier  transform  is  the  sine  function 

F[g(t)]  =  G{co)  =  ^sin^Lj  =  Tsinc^j  (A-15) 

To  obtain  a  finite  train  of  sinewave  pulses,  we  multiply  the  periodic  pulse  train  by  the 
gating  function.  Next,  we  invoke  the  frequency  convolution  theorem  to  obtain  the  Fourier 
transform  of  a  finite  train  of  sinewave  pulses 

wr(0]  =  ■  g(t )]  =  -^-F{(0)  ®  G(0))  =  —  f  F(V)G(0)-  v)dv  (A-16) 

2  7t  2  K 

Convolving  F(co)  with  the  transform  of  the  gate  function, G(ffl),  in  the  frequency  domain 
yields 
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1  I  n=+°°  I 

F[fr(r)]  =  —F(Q))®G(Q})=\^cnS(O)-nQ}0)\®G(co) 
In 


n=-H“ 


=  YjCnS((O-n(O0)®G{(O) 

n~-oo 

=  ^cnG(co-no) o) 

«=+«» 

=  '£c_nG(0)+n(0 o) 

n=+°° 

=  ^c*»G(<y+n&>0) 


(A-17) 


Substituting  Eq.  A- 15  into  A-17,  we  finally  obtain 

ET  f(tW-T  Iv”  r  ^ll 

F  /r  Lti  "  ( O)+nco0)T/2 


(A-18) 


0 ;  «=+0° 

n/,«)]=fiS 


*0  *= 


-ft 

*  n  n— _ no 


*0  n=- 


sin[(flc -n^Nn  /  fflc]sin[(<fl  +nffl0)T  /  2] 
(£yc-n6>0)(ft>+«tf>0) 

sin[(ft?c  +  wftJn)Mr/  fl7c]sin[(6J+nfi;0)r  /  2] 
Oc  +«t»0)(t»+n<2>0) 


Z 


/I=-oo 


sin[(6Jc  -n<2)0)T/2]sin[(6J+nft)0)r  /  2] 
~  (<»c-/tiW0)(ty+n670) 

sin[(<2?r  +  n67n)T/2]sin[(6?+nfi)0)r  /  2] 
(6>c  +n<»0)(ft?+n<y0) 


(A-19) 


The  Fourier  transform  of  the  periodic  sinewave  pulse  train  of  finite  duration  T  (as  are  all 
physical  signals)  is  the  sum  of  terms  given  by  Eq.  A-19,  where  <a0  is  the  repetition  rate  of 
the  pulses  in  the  train.  Each  term  of  this  sum  can  be  plotted  as  an  individual  amplitude 
distribution  centered  at  co  =  nco0 .  As  the  duration  T  increases,  the  amplitude  distributions 
become  more  compact  in  the  frequency  dimension.  As  T  goes  to  infinity,  the  total 
amplitude  distribution  approaches  a  line  spectrum. 

Thus,  four  frequency  scales  are  required  to  completely  describe  a  finite  train  of  pulsed 
sinewaves:  the  bandwidth  of  the  spectral  lines  (1  /  7);  spacing  of  the  spectral  lines 
(PRF  =  1  /  ?0 ) ;  bandwidth  of  the  sine  function  envelope  (1  / 1);  and  the  carrier  frequency 

fc  • 
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Appendix  B  Gaussian  spectrum 

In  this  appendix  we  derive  the  Fourier  transform  of  an  amplitude  modulated  cosine  carrier 
where  the  modulation  waveform  is  the  Gaussian  envelope.  Let  a  Gaussian  pulsed 
sinusoid  of  frequency  Q)0,  amplitude  A  and  pulse  width  V2?  be  represented  by  the  real 
part  of 

f(t)  =  Ae-‘2,(2x)  2ei(0°‘  .  (B-l) 

where  V2?  =  standard  deviation  .  Since  this  function  is  symmetric  with  respect  to  the 
origin  (even),  its  spectral  function  F(  0) )  is  also  even.  We  calculate  its  spectral  density 
function  as  follows 


f(t)~  A- Re 


exp 


-r 

(2t) 


-+ico0t 


-  Ae 


(2t) 


T  I  e'Vor 


(B-2) 


Taking  the  Fourier  transform  of  Eq.  B-2,  we  have 


A  °° 

F{(0)  =  —  J  exp 


-iox  +  ico^t  — 


At 


,  A  r 
dt+—  J  exp 


-iOX -i(L)J - - 

At 


dt 


(B-3) 
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Completing  the  square  by  adding  and  subtracting  [it(co  -  ©0)]  inside  the  exponent  of  the 
first  integral,  and  similarly  adding  and  subtracting  [jt(co  +  ©0)]  inside  the  exponent  of  the 
second  integral,  we  obtain 


A  r 

F(G))  =  —  J  exp 
2  J 


A  r 

--jexp 


-I  -+it(co-o)0)  I  -t2(a)-(o0y 

It  ) 

\  2 

-\—  +  it((Q  +  CO0)\  -t2(OJ+O}0)2 

2 1  J 


dt 


(B-4) 


dt 


Next,  we  make  a  change  of  variable  in  each  integral.  Designate  the  first  squared  term  in 
the  exponential  as  x,  then  dt  =  2tdx.  After  making  this  substitution  the  limits  of 
integration  remain  infinite,  and  we  have 

FW^Ate-*1^1  je^dx+Ate-^^je^dx  (B-5) 

—00  — 00 

From  the  integral  tables  we  find  that 
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So  Eq.  B-5  simplifies  to 


F((0)=4KAK-^"-“’f  +VffAnT'W*>!  (B-7) 

which  is  the  sum  of  two  Gaussian  curves  shifted  up  and  down  by  the  angular  frequency  of 
the  carrier.  If  the  original  pulse,  Eq.  B-l,  is  narrowed  by  making  T  smaller,  the  spectral 
density  curve,  Eq.  B-7,  becomes  broader  with  reduced  amplitude  exhibiting  the  inverse 
relationship  between  temporal  width  and  spectral  width  in  accordance  with  the  scaling 
property  of  the  Fourier  transform: 

f{th)t==^\i\F{TCC)  (B-8) 

Finally,  we  note  that  the  Fourier  transform  of  the  Gaussian  modulated  sinewave  is  real¬ 
valued. 


Appendix  C  FORTRAN  codes  sinewave-Debye  case 


program  SDT 

****************************************************************** 

*  SDT:  Propagates  a  finite  train  of  sinewave  pulses  into  a  Debye 

*  half-space  using  the  Fourier  transform  and  the  Fresnel  trarts- 

*  mission  coefficient  for  normal  incidence  (air/water  interface). 

* 

*  ...  Z  FIXED,  TIME  VARIES... 

* 

*  by  J.Franzen  7-12-95 


parameter(  npts  =  4000,  mpts  =  2000 ) 

complex*16  j,  FOUR(npts),  ZERO,  k,  Et,  eps,  DENOM,  ior,  E 

implicit  real*8  (a-h,o-z) 

logical  ex 

pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
ZERO  =  dcmplx(0d0,0d0) 
j  =  dcmplx(0d0,1d0) 

inquire  (file=’SD_spect.dat’,exist=ex) 
if  (ex)  then 

open(9,file=’SD_spect.dat’,status=’OLD’) 
rewind  9 
else 

open(9,file=’SD_spect.dat’,status=’NEW’) 
end  if 

fc  =  IdO  !  1  GHz  Carrier 

fmax  =  1d1  !  (+/-)  10  GHz  Integration  Limit 

tau  =  1  dl  !10  nsec  On-time  =  1 0  Cycles 

prf  =  1d0/2d1  !  PRF  =  1/PRI  =  1/20  =  50%  Duty  Cycle 

dwell  =  1d2  !  FIVE  Pulses 

*  dwell  =  1d1  !  ONE  Pulse 

df  =  fmax/dfloat(npts) 

print*, ’Repetition  period  =  ’,1d0/df,’ nanosec’ 
dw  =  twopi*df 

Check  if  Fourier  spectrum  already  exists;  Jf-so,  read  into  array  FOUR(n). 
if  (ex)  then 

print*, ’Reading  Fourier  spectrum...’ 
do  n=1,npts 
read  (9,200)  f,  FOUR(n) 
end  do 
print*, ’Done!’ 
goto  1 
end  if 

print*, ’Begin  ’,npts,’  point  Fourier  transform  now....’ 
call  fourier(npts,  fmax,  fc,  tau,  prf,  dwell,  FOUR,  dw) 
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print*, ’End  of  Fourier  transform  computation!’ 


do  n=1  ,npts 
f  =  (n-1)*df 

write(9,200)  f,  dreal(FOUR(n»,  dimag(FOUR(n)) 

200  f  ormat(1  x,f  1 0.5, 2x1  pE20.1 2,2x1  pE20.1 2) 

end  do 

close  (unit=  9) 

1  CONTINUE 

write(*,’(”  Please  enter  z  (m):  ”,$)’) 
read*,z 

write(*,’(”  Plz  enter  t-min  (ns): 
read*,tmin 

write(*,’(”  Plz  enter  t-max  (ns):  ”,$)’) 
read*,tmax 

inquire  (file=’SDT.dat’,exist=ex) 
if  (ex)  then 

open(10,fi!e=’SDT.dat’,status=’OLD’) 
rewind  10 
else 

open(10,file=’SDT.dat’,status=’NEW’) 
end  if 

c  =  .2997d0  !  light  speed  (m/ns) 

tau  =  8.1  d-3  !  relaxation  const,  scaled  to  19.648758  GHz  loss  peak 

eps_0  =  1  d7/(4d0*pi*c**2)  !  Permittivity  of  free  space 
it  =  0 

dt  =  (tmax-tmin)/dfloat(mpts) 
do  t  =  tmin,  tmax,  dt 
Et  =  ZERO 
do  n  =  1,npts 
w  =  dfloat(n-1)*dw 
DENOM  =  dcmplx(1d0,  w*tau) 
eps  =  5.5d0  +  72.7dO/DENOM 
if  (w  .ne.  OdO)  then 

eps  =  eps  -  j*(1  d-5)/w/eps_0  ! ...  subtract  ionic  conductivity 

end  if 

ior  =  cdsqrt(eps)  v  . 

*  k  =  w/c*ior 

k  =  w/c  !  •••  no  dispersion,  for  check  out 

E  =  FOUR(n)*cdexp(-j*(k*z  -  w*t)) ! ...  no  dispersion,  for  check  out 

*  E  =  FOUR(n)*2dO/(1dO+ior)*cdexp(-j*(k*z  -  w*t)) 

Et  =  Et  +  E 

*  Trapezoidal  rule  end-point  correction: 

if(n.eq.1  .or.  n.eq.npts)  then 
Et  =  Et  -  .5dO*E 
end  if 
end  do 
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it  =  it  +  1 

write(1 0,200)  t,dreal(Et)/twopi*dw 

*  write(*  ,  * )  t,dreal(Et)/twopi*dw 
if(mod(it,100).eq.0)  then 

write(*  ,*)’ it  =’, it 
end  if 
end  do 

close  (unit=10) 

stop 

END 

subroutine  fourier(  npts,  fmax,  fc,  tau,  prf,  dwell,  four,  dw  ) 
*********************************************************************** 

*  Computes  the  Fourier  transform  of  a  FINITE  train  of  pulsed  sine- 

*  waves. 

*  Modified  to  compute  spectrum  for  positive  frequencies  only. 

* 

*  Inputs: 

*  npts  -  Number  of  points  at  which  to  compute  the  Fourier  spectrum 

*  fmax  -  Maximum  frequency  range  of  Fourier  spectrum  (Hz) 

*  fc  -  Carrier  frequency  of  pulse  train  (Hz) 

*  tau  -  Pulse  on-time  (sec) 

*  prf  -  Pulse  repetition  frequency  (Hz) 

*  dwell  -  Duration  of  pulse  train  (sec) 

* 

*  Outputs: 

*  four  -  Complex  array  of  size  npts  containing  Fourier  spectrum 

*  dw  -  Frequency  resolution:  fmax/npts*twopi  (rad/sec) 

* 

*  by  J.Franzen  5-24-95 

*********************************************************************** 

IMPLICIT  REAL*8  (a-h,o-z) 

COMPLEX*16  j,  ZERO,  FOUR(npts) 

sinc(q)  =  dsin(q)/q 
pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
j  =  dcmplx(0d0,1d0) 

ZERO  =  dcmplx(0d0,0d0) 
do  n  =  1,  npts 
FOUR(n)  =  ZERO 

end  do  , 

omegaj)  =  2d0*pi*prf 
omega_c  =  2d0*pi*fc 
t_0  =  IdO/prf 
T  =  dwell 

dw  =  twopi*fmax  /  dfloat(npts) 
do  n  =  1,  npts 

w  =  dfloat(n-1)*dw 

do  L  = -1000, 1000, 1 
sum  =  w  +  L*omega_0 
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sum  =  sum  *  T/2d0 
dif  =  omega_c  -  L*omega_0 
dif  =  dif  *  tau/2d0 
if  (sum  .eq.  OdO)  then 
ap  =  tau*T/4dO*sinc(dif) 
else  if  (dif  .eq.  OdO)  then 
ap  =  tau*T/4dO*sinc(sum) 
else  if  (sum  .eq.  OdO  .and.  dif  .eq.  OdO)  then 
ap  =  tau*T/4dO 
else 

ap  =  tau*T/4dO*sinc(dif)*sinc(sum) 
end  if 

suml  =  omega_c  +  L*omega_0 
suml  =  suml  *  tau/2d0 
if  (sum  .eq.  OdO)  then 
an  =  tau*T/4dO*sinc(sum1) 
else  if  (suml  .eq.  OdO)  then 
an  =  tau*T/4dO*sinc(sum) 
else  if  (sum  .eq.  OdO  .and.  suml  .eq.  OdO)  then 
an  =  tau*T/4dO 
else 

an  =  tau*T/4dO*sinc(sum1  )*sinc(sum) 
end  if 

FOUR(n)  =  FOUR(n)  +  ap  -  an 
end  do 

FOUR(n)  =  4dO*j/t_0*FOUR(n) 
if  (mod(n,100).eq.0)  then 
print*, ’n  =  ’,n 
end  if 
end  do 

RETURN 

END 
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program  SDT2 

****************************************************************** 

*  SDT2:  Propagates  a  finite  train  of  sinewave  pulses  into  a  Debye 

*  half-space  using  the  Fourier  transform  and  the  Fresnel  trans- 

*  mission  coefficient  for  normal  incidence  (air/water  interface). 

* 

*  ...  Z  FIXED,  TIME  VARIES... 

* 

*  ...  MODIFIED  to  remove  peroidicity 

* 

*  by  J.Franzen  2-05-96 
****************************************************************** 

parameter(  npts  =  4000,  mpts  =  2000 ) 

complex*16  j,  FOUR(npts),  ZERO,  k,  Et,  eps,  DENOM,  ior,  E 

implicit  real*8  (a-h,o-z) 

logical  ex 

pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
ZERO  =  dcmplx(0d0,0d0) 
j  =  dcmplx(0d0,1d0) 

inquire  (file=’SD_spect.dat’,exist=ex) 
if  (ex)  then 

open(9,file=’SD_spect.dat’,status=’OLD’) 
rewind  9 
else 

open(9,file=’SD_spect.dat’,status=’NEW’) 
end  if 

fc  =  IdO  !  1  GHz  Carrier 

fmax  =  1d1  !  (+/-)  10  GHz  Integration  Limit 

tau  =  1d1  !  10  nsec  On-time  =  10  Cycles 

prf  =  1  d0/2d1  !  PRF  =  1/PRI  =  1/20  =  50%  Duty  Cycle 

dwell  =  1d2  !  FIVE  Pulses 

*  dwell  =  1d1  I  ONE  Pulse 

df  =  fmax/dfloat(npts) 
period  =  1  dO/df 

print*, ’Repetition  period  =  ’, period,’  nanosec’ 
dw  =  twopi*df 

Check  if  Fourier  spectrum  already  exists^lf  so,  read  into  array  FOUR(n). 
if  (ex)  then 

print*, ’Reading  Fourier  spectrum...’ 
do  n=1  ,npts 
read  (9,200)  f,  FOUR(n) 
end  do 
print*, ’Done!’ 
goto  1 
end  if 

print*, ’Begin  ’,npts,’  point  Fourier  transform  now....’ 
call  fourier(npts,  fmax,  fc,  tau,  prf,  dwell,  FOUR,  dw) 
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print*, ’End  of  Fourier  transform  computation!’ 


do  n=1  ,npts 
f  =  (n-1)*df 

write(9,200)  f,  dreal(FOUR(n)),  dimag(FOUR(n)) 

200  format(1  x,f  1 0.5, 2x1  pE20.1 2,2x1  pE20.1 2) 

end  do 

close  (unit=  9) 

1  CONTINUE 

write(*,’(”  Please  enter  z  (m):  ”,$)’) 
read*,z 

write(*,’(”  Plz  enter  t-min  (ns):  ”,$)’) 
read*,tmin 

write(*,’(”  Plz  enter  t-max  (ns):  ”,$)’) 
read*,tmax 

inquire  (file=’SDT.dat’,exist=ex) 
if  (ex)  then 

open(10,file=’SDT.dat’,status=’OLD’) 
rewind  10 
else 

open(10,file=’SDT.dat’,status=’NEW’) 
end  if 

c  =  .2997d0  !  light  speed  (m/ns) 

tau  =  8.1  d-3  !  relaxation  const,  scaled  to  19.648758  GHz  loss  peak 

eps_0  =  1  d7/(4d0*pi*c**2)  !  Permittivity  of  free  space 
it  =  0 

dt  =  (tmax-tmin)/dfloat(mpts) 
do  t  =  tmin,  tmax,  dt 
Et  =  ZERO 
do  n  =  1,npts 
w  =  dfloat(n-1)*dw 
DENOM  =  dcmplx(1d0,  w*tau) 
eps  =  5.5d0  +  72.7dO/DENOM 
if  (w  .ne.  OdO)  then 

eps  =  eps  -  j*(1  d-5)/w/eps_0  ! ...  subtract  ionic  conductivity 

end  if 

ior  =  cdsqrt(eps) 
k  =  w/c*ior 

*  k  =  w/c  !  •••  no  dispersion,  for  check  out 

*  E  =  FOUR(n)*cdexp(-j*(k*z  -  w*t)) ! ...  no  dispersion,  for  check  out 

*  E  =  FOUR(n)*2dO/(1dO+ior)*cdexp(-j*(k*z  -  w*t)) 

E  =  U(t)*FOUR(n)*2dO/(1dO+ior)*cdexp(-j*(k*z  -  w*t))  - 
&U(t-period+dwell/2.)*FOUR(n)*2dO/(1dO+ior)* 
&cdexp(-j*(k*z-w*(t-period))) 

Et  =  Et  +  E 

*  Trapezoidal  rule  end-point  correction: 
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if(n.eq.1  .or.  n.eq.npts)  then 
Et  =  Et  -  .5dO*E 
end  if 
end  do 

it  =  it  +  1 

write(1 0,200)  t,dreal(Et)/twopi*dw 

*  write(*  ,  * )  t,dreal(Et)/twopi*dw 
if(mod(it,100).eq.0)  then 

write(*,*)  ’  it  =  %it 
end  if 
end  do 

close  (unit=10) 

stop 

END 

subroutine  fourier(  npts,  fmax,  fc,  tau,  prf,  dwell,  four,  dw ) 
*********************************************************************** 

*  Computes  the  Fourier  transform  of  a  FINITE  train  of  pulsed  sine- 

*  waves. 

*  Modified  to  compute  spectrum  for  positive  frequencies  only. 

* 

*  Inputs: 

*  npts  -  Number  of  points  at  which  to  compute  the  Fourier  spectrum 

*  fmax  -  Maximum  frequency  range  of  Fourier  spectrum  (Hz) 

*  fc  -  Carrier  frequency  of  pulse  train  (Hz) 

*  tau  -  Pulse  on-time  (sec) 

*  prf  -  Pulse  repetition  frequency  (Hz) 

*  dwell  -  Duration  of  pulse  train  (sec) 

* 

*  Outputs: 

*  four  -  Complex  array  of  size  npts  containing  Fourier  spectrum 

*  dw  -  Frequency  resolution:  fmax/npts*twopi  (rad/sec) 

* 

*  by  J.Franzen  5-24-95 

*********************************************************************** 

IMPLICIT  REAL*8  (a-h,o-z) 

COMPLEX*16  j,  ZERO,  FOUR(npts) 

sinc(q)  =  dsin(q)/q 
pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
j  =  dcmplx(0d0,1d0) 

ZERO  =  dcmplx(0d0,0d0) 
do  n  =  1 ,  npts 
FOUR(n)  =  ZERO 
end  do 

omega_0  =  2d0*pi*prf 
omega_c  =  2d0*pi*fc 
t_0  =  IdO/prf 
T  =  dwell 

dw  =  twopi*fmax  /  dfloat(npts) 
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do  n  =  1,  npts 

w  =  dfloat(n-1)*dw 

do  L  = -1000, 1000, 1 
sum  =  w  +  L*omega_0 
sum  =  sum  *  T/2d0 
dif  =  omega_c  -  L*omega_0 
dif  =  dif  *  tau/2d0 
if  (sum  .eq.  OdO)  then 
ap  =  tau*T/4dO*sinc(dif) 
else  if  (dif  .eq.  OdO)  then 
ap  =  tau*T/4dO*sinc(sum) 
else  if  (sum  .eq.  OdO  .and.  dif  .eq.  OdO)  then 
ap  =  tau*T/4dO 
else 

ap  =  tau*T/4dO*sinc(dif)*sinc(sum) 
end  if 

suml  =  omega_c  +  L*omega_0 
suml  =  suml  *  tau/2d0 
if  (sum  .eq.  OdO)  then 
an  =  tau*T/4dO*sinc(sum1) 
else  if  (suml  .eq.  0d0)  then 
an  =  tau*T/4dO*sinc(sum) 
else  if  (sum  .eq.  OdO  .and.  suml  .eq.  OdO)  then 
an  =  tau*T/4dO 
else 

an  =  tau*T/4dO*sinc(sum1  )*sinc(sum) 
end  if 

FOUR(n)  =  FOUR(n)  +  ap  -  an 
end  do 

FOUR(n)  =  4dO*j/t_0*FOUR(n) 
if  (mod(n,100).eq.0)  then 
print*, ’n  =  ’,n 
end  if 
end  do 

RETURN 

END 

real*8  function  U(t) 

if  (t.ge.OdO)  then  >.  . 

U  =  IdO 
else 
U  =  OdO 
end  if 
return 
end 
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Appendix  D  FORTRAN  code  Gaussian-Lorentz  case. 


program  GLT2 

********************************************************************* 

*  GLT2:  Propagates  a  single  Gaussian  pulse  into  a  Lorentzian  half- 

*  space  using  the  Fourier  transform  and  the  Fresnel  transmission 

*  coefficient  for  normal  incidence. 

* 

****  Modified  to  use  Oughstun’s  nondimensional  space-time  parameter 
****  THETA 

********************************************************************* 

* 

*  ...  Z  FIXED,  TIME  VARIES  ... 

* 

*  by  J.Franzen  9-29-95 
********************************************************************* 

parameter(  npts  =  2**13,  mpts  =  10  000 ) 
implicit  real*8  (a-h,o-z) 

complex*16  j,  FOUR(npts),  ZERO,  Et,E,  eps,  DENOM 

complex*16T12(npts),k(npts),ior(npts) 

logical  ex 

pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
ZERO  =  dcmplx(OdO.OdO) 
j  =  dcmplx(0d0,1d0) 

c  =  2.997d8  /  1d16  !  light  speed  (m/deci-femto-sec) 

*  Set  values  for  demonstration: 

f c  =  .91 5d0  !  Carrier  Freq;  X  10**16  Hz 

fmax  =  1d1  !  Integration  Limit;  X  10**16  Hz 

tau  =  .25d0  !  On-time;  X  10**(-1 6)  sec 

df  =  fmax/dfloat(npts) 

print*, ’Repetition  period  =  ’,1  d-1/df ,’  femtosec’ 
dw  =  twopi*df 

Check  if  Fourier  spectrum  already  exists.  If  so,  read  into  array  FOUR(n). 
inquire  (file=’GLT2_spect.dat’,exist=ex) 
if  (ex)  then 

open(9,file=’GLT2_spect.dat’,status=’OLD’) 
rewind  9  v- 

else 

open(9,file=’GLT2_spect.dat’,status=’NEW’) 
end  if 

if  (ex)  then 

print*, ’Reading  Fourier  spectrum...’ 
do  n=1  ,npts 
read  (9,200)  f,  FOUR(n) 
end  do 
print*, ’Done!’ 
goto  1 
end  if 
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print*, ’Begin  ’,npts,’  point  Fourier  transform  now....’ 

call  gfour(npts,  fmax,  fc,  tau,  FOUR,  dw) 

print*, ’End  of  Fourier  transform  computation!’ 

do  n=1,npts 
f  =  (n-1)*df 

write(9,200)  f,  FOUR(n) 

200  format(1  x,f  1 0.5, 2x1  pE20.1 2,2x1  pE20.1 2) 

end  do 

close  (unit=  9) 

1  CONTINUE 

Compute  Refl.,  Trans,  coefs  &  k-numbers  for  LORENTZ  half-space, 
omegaj)  =  4d0  !  Resonant  Frequency;  X  10**16  rad/s 

delta  =  .28d0  !  Damping  Constant;  X  10**16  rad/s 

do  n=1  ,npts 

w  =  dfloat(n-1)*dw 

DENOM  =  w**2  -  2d0*j*w*delta  -  omega_0**2 
eps  =  IdO  -  (1.25dO*omega_0**2)/DENOM 
ior(n)  =  cdsqrt(eps) 

T1 2(n)  =  2d0/(  1  dO  +  ior(n) )  !  TE  Polarization 

k(n)  =  w/c*ior(n) 
end  do 

write(*,’(”  Please  enter  z  (m):  ”,$)’) 
read*,z 

print*, ’z  (microns)  =  ’,z*1d6 

write(*, ’(”  Plz  enter  t-min  (fs):  ”,$)’) 

read*,tmin 

print*, ’t-min  =  ’,tmin 

tmin  =  tmin*1d1 

write(*,’(”  Plz  enter  t-max  (fs):  ”,$)’) 
read*,tmax 

print*, ’t-max  (fs)  =  ’,tmax 
tmax  =  tmax*1d1 

inquire  (file=’GLT2.dat’,exist=ex) 

if  (ex)  then  >.  . 

open(10,file=’GLT2.dat’,status=’OLD’) 
rewind  10 

else 

open(10,file=’GLT2.dat’,status=’NEW’) 
end  if 

*  For  each  value  of  z,  compute  the  inverse  Fourier  integral  by 

*  numerical  quadrature  using  the  trapezoidal  rule. 

it  =  0 

dt  =  (tmax-tmin)/dfloat(mpts) 
do  t  =  tmin,  tmax,  dt 
theta  =  c*t/z 
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it  =  it  +  1 
Et  =  ZERO 
do  n  =  1  ,npts 

w  =  dfloat(n-1)*dw 

*  E  =  FOUR(n)*T12(n)*cdexp(-j*(k(n)*z  -  w*t) ) 

E  =  FOUR(n)*T12(n)*cdexp(-j*z/c*w*(ior(n)-theta)) 

Et  =  Et  +  E  !  Transmitted  Field  only 

*  Trapezoidal  rule  end-point  corrections: 

if(n.eq.1  .or.  n.eq.npts)  then 
Et  =  Et  -  ,5dO*E 
end  if 


end  do  I  End  of  w-loop.... 

Et  =  Et/twopi*dw 

if(mod(it,100).eq.0)  then 
write(Y)  ’  it  =  \it 
end  if 

write(1 0,201 )  theta,  -dreal(Et) 

201  format(1  x,2x1  pE20.1 2,2x1  pE20.1 2) 

end  do  !  End  of  t-loop.... 

close  (unit=10) 

stop 

END 

subroutine  gfour(  npts,  fmax,  fc,  tau,  *ou^^2»******* 

*  Computes  the  Fourier  transform  of  a  Gaussian  pulse. 

* 

*  Inputs: 

*  npts  -  Number  of  points  at  which  to  compute  Fourier  spectrum 

*  fmax  -  Maximum  frequency  range  of  Fourier  spectrum  (Hz) 

*  tau  -  sigma/sqrt(2) 

* 

*  Outputs: 

*  four  -  Complex  array  of  size  npts  containing  Fourier  spectrum 

*  dw  -  Frequency  resolution:  fmax/npts*twopi  (rad/sec) 


IMPLICIT  REAL*8  (a-h,o-z) 
COMPLEX*16  j,  ZERO,  FOUR(npts) 

pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
j  =  dcmplx(0d0,1d0) 

ZERO  =  dcmplx(0d0,0d0) 

omega_c  =  2d0*pi*fc 

dw  =  twopi*fmax  /  dfloat(npts) 

do  n  =  1,  npts 
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w  =  dfloat(n-1)*dw 

FOUR(n)  =  dexp(-(tau**2)*(w-omega_c)**2) 
FOUR(n)  =  j*dsqrt(2dO*twopi)*tau*FOUR(n) 
print*, ’four(’,n,’)  =  ’,four(n) 
if  (mod(n,100).eq.0)  then 
print*, ’n  =  ’,n 
end  if 
end  do 


RETURN 


Appendix  E  FORTRAN  code  water  impulse  response 


program  IMP22 


*  IMP22:  Performs  time-domain  convolution  on  sinewave  pulse  using 

*  the  impulse  response  of  water. 

* 


*  by  J.Franzen  2-01-96 

****************************************************************** 

parameter(  npts  =  4000,  mpts  =  8000 ) 
complex*16  j,  ZERO,  k,  eps,  DENOM,  ior,  E 
implicit  real*8  (a-h,o-z) 

real*8  IMP(mpts),tau(mpts),X(npts),Y(npts+mpts-1) 
logical  ex 

pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
ZERO  =  dcmplx(0d0,0d0) 
j  =  dcmplx(0d0,1d0) 

fmax  =  2d9  !  (+/-)  20  GHz  Integration  Limit 

df  =  fmax/dfloat(npts) 
dw  =  twopi*df 

write(Y(”  Please  enter  z  (m):  ”,$)’) 
read*,z 

Check  if  Impulse  response  already  exists.  If  so,  read  into  array  IMP(m). 
inquire  (file=’IMP22.dat’,exist=ex) 
if  (ex)  then 

open(10,file=’IMP22.dat’,status=’OLD’) 
rewind  10 
else 

open(10,file=’IMP22.dat’,status=’NEW’) 
end  if 

if  (ex)  then 

print*, ’Reading  impulse  response...’^ 
do  m=1  ,mpts 

read  (10,200)  tau(m),IMP(m) 
end  do 
print*, ’Done!’ 
goto  1 
end  if 

Compute  Impulse  Response ... 

c  =  .2997d9  !  Light  speed  (m/s) 

taw  =  8.1  d-1 2  I  Water  relaxation  const. 
eps_0  =  8.854d-12  I  Permittivity  of  free  space 
it  =  0 
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dT  =  1dO/(2dO*fmax) 
do  m  =  1,mpts 
IMP(m)  =  OdO 
do  n  =  1  ,npts 
w  =  dfloat(n-1)*dw 
DENOM  =  dcmplx(1d0,  w*taw) 
eps  =  5.5d0  +  72.7dO/DENOM 
if  (w  .ne.  OdO)  then 

eps  =  eps  -  j*(1d-5)/w/eps_0  ! ...  subtract  ionic  conductivity 

end  if 

ior  =  cdsqrt(eps) 
k  =  w/c*ior 

*  k  =  w/c  ! ...  no  dispersion,  for  check  out 

*  E  =  cdexp(-j*(k*z  -  m*w*dT))  ! ...  no  dispersion,  for  check  out 

E  =  2d0/(1  dO+ior)*cdexp(-j*(k*z  -  m*w*dT)) 

*  IMP(m)=IMP(m)+dreal(E)*dcos(m*pi*dT)-dimag(E)*dsin(m*pi*dT) 

IMP(m)  =  IMP(m)  +  dreal(E) 
end  do 

it  =  it  + 1 
tau(m)  =  m*dT 

IMP(m)  =  2dO*IMP(m)/twopi*dw 
write(1 0,200)  tau(m),IMP(m) 

200  f  ormat(1  x,1  pel  4.6,2x,1  pel  4.6) 
if(mod(it,100).eq.0)  then 
write(Y)  ’  it  =  ’,it 
end  if 
end  do 

close  (unit=10) 

1  CONTINUE 

Compute  the  Convolution  Integral: 

inquire  (file=’PULSE22.dat’,exist=ex) 
if  (ex)  then 

open(1 1  ,file=’PU  LSE22.dat’,  status=’OLD’) 
rewind  11 
else 

open(11,file=’PULSE22.dat’,status=’NEW’) 
end  if 

fc  =  1d9  !!!  1  GHz  Carrier  .Frequency 

on_time  =  Id-8 !!!  10  ns  on-time  (10 cycles) 

it  =  0 

dT  =  1d0/(2d0*fmax) 
do  n=1  ,npts 
t  =  n*dT 

X(n)  =  pulse(t,  fc,on_time) 
end  do 

call  CONVL  (X,  npts,  IMP,  mpts,  Y,  ly) 
do  n=1  ,ly 

write(1 1,200)  n*dT,  Y(n)*dT 
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end  do 

close  (unit=11) 

stop 
END 

real*8  function  pulse  (t,  fc,on_time) 
implicit  real*8  (a-h,o-z) 
pi  =  dacos(-ldO) 

if  (t  .It.  on_time)  then 
pulse  =  dsin(2d0*pi*fc*t) 
else 

pulse  =  OdO 
end  if 

return 
END 

subroutine  CONVL  (X,  n,  H,  m,  Y,  ly) 

Subroutine  CON VL  computes  the  convolution 
between  the  sequences  X  and  H 
Parameters 

X  :  array  containing  sequence  X 
n  :  length  of  sequence  X 
H  :  array  containing  sequence  H 
m  :  length  of  sequence  H 

Y  :  array  containing  the  convolution  of  X  and  H 
ly  :  length  of  the  sequence  Y 

real*8  X,H,Y 

dimension  X(1),H(1),Y(1) 

ly  =  n  +  m  - 1 
do  k  =  1  ,  ly 
Y(k)  =  OdO 
end  do 

do  1  i  =  1  ,  n 
do  1 j  =  1  , m 

k  =  i  +  j  - 1  v  ,.r. 

1  Y(k)  =  Y(k)  +  X(i)  *  H(j) 

return 

END 
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Appendix  F  FORTRAN  code  Lorentz  medium  impulse  response 


program  IMP3 

********************************************************************* 

*  IMP3:  Performs  time-domain  convolution  on  sinewave  pulse  using  the 

*  impulse  response  of  a  single  absorption  band  Lorentz  modim. 

* 

*  1  Pulse  @  8d15  Hz  carrier,  10  cycles-on 


*  by  J.Franzen  10-23-95 
********************************************************************* 

parameter  npts  =  2**13,  mpts  =  2**14 ) 
implicit  real*8  (a-h,o-z) 
complex*16  j,  ZERO,  E,  eps,  DENOM 
complex*16  T12(npts),k(npts),ior 
real*8  IMP(mpts),tau(mpts),X(npts),Y(npts) 
logical  ex 

pi  =  dacos(-ldO) 
twopi  =  2d0*pi 
ZERO  =  dcmplx(0d0,0d0) 
j  =  dcmplx(0d0,1d0) 

c  =  2.997d8  / 1  dl  6  !  light  speed  (m/deci-f emto-sec) 

*  Set  values  for  demonstration: 

fc  =  .8d0  !  Carrier  Freq;  X  10**16  Hz 

fmax  =  2d1  !  Integration  Limit;  X  10**16  Hz 

taw  =  12.5d0  !  On-time;  X  1 0**(-1 6)  sec 

df  =  fmax/dfloat(npts) 

print*, ’Repetition  period  =  \1d-1/df,’  femtosec’ 
dw  =  twopi*df 

Check  if  impulse  response  already  exists.  If  so,  read  into  array  IMP(m). 
inquire  (file=’IMP3.dat’,exist=ex) 
if  (ex)  then 

open(9,file=’IMP3.dat’,status=’OLD’) 
rewind  9 
dSG 

open(9,file=’IMP3.dat’,status=’NEW’) 
end  if 

if  (ex)  then 

print*, ’Reading  impulse  rsponse 
do  n=1  ,mpts 

read  (9,200)  tau(m),IMP(m) 
end  do 
print*, ’Done!’ 
goto  1 
end  if 
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Compute  impulse  response ... 

omega_0  =  4d0  !  Resonant  Frequency;  X  10**16  rad/s 

delta  =  .28d0  !  Damping  Constant;  X  10**16  rad/s 

do  n=1,npts 

w  =  dfloat(n-1)*dw 

DENOM  =  w**2  -  2d0*j*w*delta  -  omega_0**2 
eps  =  IdO  -  (1.25dO*omega_0**2)/DENOM 
ior  =  cdsqrt(eps) 

T12(n)  =  2d0/(  IdO  +  ior ) 
k(n)  =  w/c*ior 
end  do 

write(*,’(”  Please  enter  z  (microns):  ”,$)’) 
read*,z 

print*, ’z  (microns)  =  ’,z 
z  =  z*1d-6 


it  =  0 

dT  =  1d0/(2d0*fmax) 
do  m  =  1,  mpts 
it  =  it  +  1 
IMP(m)  =  OdO 
tau(m)  =  m*dT 
do  n  =  1,npts 

w  =  dfloat(n-1)*dw 

E  =  T1 2(n)*cdexp(-j*(k(n)*z  -  w*m*dT) ) 
IMP(m)  =  IMP(m)  +  dreal(E) 
end  do 

IMP(m)  =  2dO*IMP(m)/twopi*dw 
write(9,200)  tau(m),  IMP(m) 
if  (mod(it,1 00).eq.0)  then 
write(*,*)  ’  it  =  ’,it 
end  if 

200  f  ormat(1  x,1  pEI  4.6, 2x1  pEI  4.6) 
end  do 
close  (unit=9) 

1  CONTINUE 

Compute  the  Convolution  Integral: 

inquire  (file=’PULSE3.dat’,exist=ex) 
if  (ex)  then 

open(1 1  ,file=’PU  LSE3.dat’, status=’OLD’) 
rewind  11 
else 

open(1 1  ,file=’PU  LSE3.dat’, status=’NEW’) 
end  if 

it  =  0 

dT  =  1d0/(2d0*fmax) 
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t  =  n*dT 

X(n)  =  pulse(t,  fc,taw) 
end  do 

call  CONVL  (X,  npts,  IMP,  mpts,  Y,  ly) 


do  n=1  ,ly 

write(1 1,200)  n*dT,  Y(n)*dT 
end  do 

close  (unit=11) 

stop 

END 

real*8  function  pulse  (t,  fc,on_time) 
implicit  real*8  (a-h,o-z) 
pi  =  dacos(-ldO) 

if  (t  .It.  on_time)  then 
pulse  =  dsin(2d0*pi*fc*t) 
else 

pulse  =  OdO 
end  if 

return 

END 

subroutine  CONVL  (X,  n,  H,  m,  Y,  ly) 

Subroutine  CONVL  computes  the  convolution 
between  the  sequences  X  and  H 
Parameters 

X  :  array  containing  sequence  X 
n  :  length  of  sequence  X 
H  :  array  containing  sequence  H 
m  :  length  of  sequence  H 
Y  :  array  containing  the  convolution  of  X  and  H 
ly  :  length  of  the  sequence  Y 

real*8  X,H,Y 

dimension  X(1),H(1),Y(1)  >.  „ 

ly  =  n  +  m  - 1 
do  k  =  1  ,  ly 
Y(k)  =  OdO 
end  do 

do  1  i  =  1  ,  n 
do  1 j  =  1 , m 
k  =  i  +  j  - 1 

Y(k)  =  Y(k)  +  X(i)  *  H(j) 


return 


END 


Single  Pulse  Time-history 

Unit  Amplitude,  Pulse  Width  =  10  ns,  Carrier  Freq.  =  1  GHz 


(iu/a)  P|S!d-3 


Figure  1.  Sinewave  pulse  with  1  GHz  carrier,  10  cycles 


Figure  2.  Spectrum  of  sinewave  pulse  with  1  GHz  carrier  ,  10  cycles 
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Figure  3.  RF  pulse  in  water  at  3A  meter  depth 


(UJ/A)  Piey-3 


Figure  4.  Sinewave  pulse  train  in  water 


Five  Pulse  Burst  Spectrum 

Unit  Amplitude,  Pulse  Width  =  10  ns,  Carrier  =  1  GHz,  Duty  Factor  =  50% 


o 


Figure  5.  Spectrum  of  five  sinewave  pulses:  duty  factor  =  Vi 


Single  Pulse  into  Lorentz  Half-space 

Fourier  Integral  Method:  1, 10,  20  microns 


(iu/a)  Piey-3 


(ui/a)  P|3!d-3 


(uj/a)  P|S!d-3 


Figure  7.  Single  pulse  into  Lorentzian  half-space 


Brillouin  Precursor:  Distance  vs  Time 
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Figure  8.  Brillouin  precursor:  distance  vs.  time 


Sommerfeld  Precursor:  Distance  vs  Time 

omeaa  0=4E+16/s,  b*b=20E+32 /s*s,  delta=0.28E+1 6 /s 
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Figure  9.  Sommerfeld  precursor:  distance  vs.  time 


Gaussian  Pulse  Spectrum 

Amplitude,  tau  =  .05  fs,  Carrier  =  5.75E+16  rad/s 


Figure  10.  Spectrum  of  single  Gaussian  pulse 


Gaussian  Pulse  into  Lorentz  Medium 

5.75E+16  rad/s,  tau  =  T/2  =  .05  fs,  Depth  =  1  micron 


Figure  11.  Gaussian  pulse  into  Lorentzian  medium 
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Microwave  Pul 
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Figure  13.  Sinewave  pulse  into  water  via  time-domain  convolution 


Impulse  Response  of  Lorentz  Medium 

Depth  =  1  Micron 


response  of  Lorentzian  medium 


Ultraviolet  Pulse  into  Lorentz  Medium 

8E15  Hz,  1.25  fs,  1  micron,  Time-domain  convolution 


(w/A)  P|9!d-3 


Figure  15.  UV  pulse  into  Lorentzian  medium  via  time-domain  convolution 
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Figure  16.  Four  cycle  pulse  reconstructed  using  10  rad/s  integration  bandwidth 
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Figure  17.  Four  cycle  pulse  reconstructed  using  5  rad/s  integration  bandwidth 
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Figure  18.  Four  cycle  pulse  reconstructed  using  2.5  rad/s  integration  bandwidth 


Five  Pulse  Time-history 

Unit  Amplitude,  Pulse  Width  =  10  ns,  Duty  Cycle  =  50%,  Carrier  =  1  GHz 


# 


(ui/A)  P|S!d-3 


Figure  Al.  Five  pulse  wave  train 


